

#include "vector"
#include "iostream"
#include "complex"
#include <time.h>

using namespace std;

typedef complex<float> Complex;

class LatticeZ {
public:
	vector<int> dim;
	unsigned long volume(int mu=-1) {
		long s=1;
		for(int k=dim.size()-1; k>mu; k--) s*=dim[k];
		return s;
	}
};


class SiteZ {
public:
	LatticeZ *ptr_lattice;
	unsigned long idx;
	SiteZ(LatticeZ &lattice, long idx=0) {
		ptr_lattice = &lattice;
		this->idx = idx;
	}
	void set(int x0=0, int x1=0, int x2=0, int x3=0, int x4=0, 
		int x5=0, int x6=0, int x7=0, int x8=0, int x9=0) {
			int n = ptr_lattice->dim.size();
			idx = 0;
			if(n>0) idx += ptr_lattice->volume(0)*x0;
			if(n>1) idx += ptr_lattice->volume(1)*x1;
			if(n>2) idx += ptr_lattice->volume(2)*x2;
			if(n>3) idx += ptr_lattice->volume(3)*x3;
			if(n>4) idx += ptr_lattice->volume(4)*x4;
			if(n>5) idx += ptr_lattice->volume(5)*x5;
			if(n>6) idx += ptr_lattice->volume(6)*x6;
			if(n>7) idx += ptr_lattice->volume(7)*x7;
			if(n>8) idx += ptr_lattice->volume(8)*x8;
			if(n>9) idx += ptr_lattice->volume(9)*x9;
	}
	int operator()(int mu) {
		return (idx / ptr_lattice->volume(mu)) % ptr_lattice->dim[mu];
	}
};

SiteZ operator+(const SiteZ &x, unsigned int mu) {
	// shift right on the lattice of one step in direction mu=+0,+1,+2,...
	mu = mu % x.ptr_lattice->dim.size();
	long idx = x.idx;
	long v = x.ptr_lattice->volume(mu);
	int n = x.ptr_lattice->dim[mu];
	long projection = (idx / v) % n;
	idx += ((projection+1) % x.ptr_lattice->dim[mu] - projection) * v;
	SiteZ y(*x.ptr_lattice,idx);  
	return y;
}

SiteZ operator-(const SiteZ &x, unsigned int mu) {
	// shift left on the lattice of one step in direction mu=-0,-1,-2,...
	mu = mu % x.ptr_lattice->dim.size();
	long idx = x.idx;
	long v = x.ptr_lattice->volume(mu);
	int n = x.ptr_lattice->dim[mu];
	long projection = (idx / v) % n;
	idx += ((projection+n-1) % n - projection) * v;
	SiteZ y(*x.ptr_lattice,idx);  
	return y;
}

class GaugeZ {
public:
	LatticeZ *ptr_lattice;
	int nc;
	vector<Complex> m;
	GaugeZ(LatticeZ &lattice, int nc) {
		ptr_lattice = &lattice;
		this->nc=nc;
		m.resize(lattice.volume()*lattice.dim.size()*nc*nc);
	}
	int ndim() const {
		return ptr_lattice->dim.size();
	}
	const Complex get(const SiteZ &x, int mu, int i, int j) const {
		int n = ndim(); 
		if(mu>=0) return m[((x.idx*n + (mu % n))*nc+i)*nc+j];
		else return conj(m[((x.idx*n + (-mu % n))*nc+j)*nc+i]);
	}
	Complex &operator()(const SiteZ &x, int mu, int i, int j) {
		int n = ndim(); 
		if(mu<0) throw "error";
		return m[((x.idx*ndim() + (mu % n))*nc+i)*nc+j];
	}
};

class PathProductZ {
public:
	LatticeZ *ptr_lattice;
	int nc;
	vector<Complex> m;
	PathProductZ(LatticeZ &lattice,int nc) {
		ptr_lattice = &lattice;
		this->nc=nc;
		m.resize(lattice.volume()*nc*nc);
	}
	const Complex get(const SiteZ &x, int i, int j) const {
		return m[(x.idx*nc+i)*nc+j];
	}
	Complex &operator()(const SiteZ &x, int i, int j) {
		return m[(x.idx*nc+i)*nc+j];
	}
};

typedef vector<int> Path;

int orgMain() {
	// input nc
	int nc=3;
	// input LatticeZ size
	LatticeZ spacetime;  
	spacetime.dim.resize(4);
	spacetime.dim[0]=spacetime.dim[1]=spacetime.dim[2]=spacetime.dim[3]=4;
	SiteZ x(spacetime);
	SiteZ y(spacetime);

	// test
	x.set(2,2,2,2);
	//cout << x(0) << x(1) << x(2) << x(3) << endl;
	x = x+2;
	//cout << x(0) << x(1) << x(2) << x(3) << endl;
	x=x-0;
	//cout << x(0) << x(1) << x(2) << x(3) << endl;
	// end test

	// input U
	// defines ndim x nc x nc matrices at each point of spacetime
	GaugeZ U(spacetime,nc); 
	for(x.idx=0; x.idx<spacetime.volume(); x.idx++)
		for(unsigned int mu=0; mu<spacetime.dim.size(); mu++)
			for(int i=0; i<nc; i++)
				for(int j=0; j<nc; j++)
					U(x,mu,i,j)=(i==j)?1.0:0.0;

	// input path
	Path path(6);
	path[0]=+1; path[1]=+1; path[2]=+2; path[3]=-1; path[4]=-1; path[5]=-2;

	// at every point of lattice compute P
	PathProductZ P(spacetime,nc);  
	Complex tmp[3][3];
	for(x.idx=0; x.idx<spacetime.volume(); x.idx++) {
		for(int i=0; i<nc; i++)
			for(int j=0; j<nc; j++)
				P(x,i,j)=U.get(x,path[0],i,j);
		y=x;
		for(unsigned int p=1; p<path.size(); p++) {
			y = y + path[p-1];
			int check = 0 ;
			for(int i=0; i<nc; i++)
				for(int j=0; j<nc; j++) {
					tmp[i][j]=0;
					for(int k=0; k<nc; k++)
						tmp[i][j]+=P.get(x,i,k)*U.get(y,path[p],k,j);
					P(x,i,j) = tmp[i][j];
					check += (tmp[i][j]==Complex(0,0))?0:1;
				}
				if(check!=nc) cout << "error " << check << endl;
		}
	}

	/*for (unsigned int i = 0; i < P.m.size(); i++ ) {
		cout << P.m[i].real() << " ";
	}*/

	// save P
	//cout << "if gets here and prints no error, it is ok\n";
	return 0;
}

// NOTES: the loops on x.idx need to be parallelized (the looping order is not relevant).
